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Abstract 

A class of quasilinear singularly perturbed boundary value problems 
with a turning point of attractive type is considered. The problems 
are solved numerically by a finite-difference scheme on a special dis¬ 
cretization mesh which is dense near the turning point. The scheme 
is a combination of the standard central and midpoint schemes and 
is practically second-order accurate. Pointwise accuracy is uniform 
in the perturbation parameter, e, and, moreover, L * 1 errors decrease 
when £ —> 0. This is achieved by the use of meshes which general¬ 
ize the piecewise equidistant Shishkin mesh. Two particular types of 
meshes are considered and compared. 

Keywords: singular perturbation, boundary value problem, finite differ¬ 
ences, Shishkin mesh. 


1 Introduction 

The following singularly perturbed boundary value problem was considered 
in Vulanovic and Lin [1]: 

— eu" — xb(x, u)v! + c(x, u) = 0, x € [—1,1], u(±l) = U±, (1) 

where £ is a small positive parameter, b > 0 and c = 0(\x\ + e) are suffi¬ 
ciently smooth functions, and U± are two given numbers. Some additional 

*This manuscript was prepared in February, 2000, but has remained unpublished and 

I decided to make it accessible through arXiv. 
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conditions were assumed as well, but we shall not list them here. It was 
shown that any solution of problem (1) has an interior layer of exponential 
type at the turning point x = 0. An appropriate numerical method was 
proposed, based on finite differences on a special discretization mesh dense 
in the layer. The error of this method was estimated in a discrete L 1 norm 
by 

M{Ve + e- N )N-\ ( 2 ) 

(In (2) and throughout the present paper, N denotes the number of mesh 
steps and M stands for any, in the sense of 0(1), positive constant indepen¬ 
dent of £ and N.) 

In the present paper, we are interested in improving the result from [1], 
This is achieved by applying a different kind of discretization mesh and a 
higher-order scheme. To simplify the presentation, we shall only consider 
a special case of problem (1), viz. the case b = b(u ) and c = 0. However, 
the same theoretical results hold for the general problem if the method and 
some conditions are modified appropriately. 

Singularly perturbed boundary value problems arise in various appli¬ 
cations, see Chang and Howes [2] for instance. Two recent books discuss 
numerical methods for these problems: Miller et al. [3] and Roos et al. [4]. 
Numerical methods for different types of turning point problems have at¬ 
tracted a considerable attention, let us only mention the papers by Berger 
et al. [5], Lin [6], Vulanovic and Farrell [7], Clavero and Lisbona [8], and Sun 
and Stynes [9]. Prior to [1], some weaker versions of problem (1) were con¬ 
sidered in Vulanovic [10], [11], and [12], where special discretization meshes 
were also used. The meshes belong to the class of explicitly constructed 
meshes, which means that they are formed before the discrete problem is 
solved. This requires sharp derivative estimates of the continuous solution. 
The same approach will be applied here. 

In general, the explicitly constructed meshes can be divided in two main 
classes, meshes of Bakhvalov type (B meshes) and meshes of Shishkin type 
(S meshes). The former were introduced in Bakhvalov [13] and later on 
generalized and simplified in Vulanovic [14], Many different types of sin¬ 
gular perturbation problems have been successfully solved numerically on 
B meshes. A B mesh is constructed by a smooth mesh generating func¬ 
tion which maps equidistant points into mesh points which are dense in the 
layer(s). In each layer, the mesh generating function corresponds to the in¬ 
verse of the function describing how the solution behaves there. On the other 
hand, S meshes (Shishkin [15]) are piecewise equidistant, and therefore much 
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simpler. An S mesh for a problem like (1) would typically consist of three 
equidistant parts: a fine part around x = 0 and two coarse parts of the mesh 
outside the layer, with the transition points between these parts located at 
±ay/elnN, where a > 0 is independent of e and of N. However, S meshes 
produce less accurate numerical results (see Vulanovic [16] for a comparison 
of B and S meshes). There are two possibilities for improving numerical 
results obtained on S meshes. One of them is to use more accurate schemes 
that may be more complicated but still easier to analyze on S meshes than 
on B meshes, and another one is to improve the mesh itself. This paper uses 
both approaches: a scheme which is practically second-order accurate and 
a modified and improved S mesh. 

Different improvements of S meshes have been considered so far. That 
by Linfi [17] uses a modification which makes the mesh more similar to a 
B mesh. The resulting mesh is not piecewise equidistant any more. The 
approach in Vulanovic [18] improves the S mesh while keeping it piecewise 
equidistant (the mesh has more equidistant parts which are constructed in 
a special way). A recent paper by Roos and Linfi [19] (see Linfi et al. [20] as 
well) provides for a unified theory which covers both S and B meshes (the 
latter only slightly different from those in [14]) and their generalizations. 
All these papers deal with non-turning point problems. 

We are not going to consider all the different types of meshes here. We 
shall only analyze the mesh from [18] and its slight modification. In this 
way, we show that meshes of this type can be applied also to turning point 
problems. The modification of the mesh from [18] consists of replacing N in 
the transition point formulas by 1/e. Such transition points in S meshes were 
briefly discussed in [19] in the non-turning point case. They are closer to 
the points marking the beginning and the end of the layer and they improve 
both theoretical and numerical results for problems of type (1). The two 
kinds of meshes give quite satisfactory results and the discretization scheme 
is easier to discuss on these piecewise equidistant meshes. This is why we are 
not going to consider B meshes here. Besides, the particular B mesh applied 
to (1) in [1] is even more complicated than some other, more standard, B 
meshes. 

The scheme which we shall use in this paper belongs to the class of hybrid 
(or switching) finite-difference schemes. The standard central scheme is used 
inside the layer, where the mesh is fine and the scheme is unconditionally 
stable, whereas a midpoint upwind scheme is used outside the layer, on the 
coarse part of the mesh. For such schemes, see Vulanovic [21] and [22], and 
more recently, Stynes and Ross [23], and Linfi [24], Of these papers, only [22] 
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and [24] deal with hybrid schemes for quasilinear problems. The scheme in 
[24] is less general, since it is constructed for a non-turning point problem. 
However, it is simpler and that is why we are going to use a very similar 
approach here. 

The paper is organized as follows. Precise assumptions on the continuous 
problem and properties of its solution are given in section 2. In addition to 
the interior turning point case, a boundary turning point is also considered. 
In that case, the turning point is still x = 0, but the interval [0,1] is consid¬ 
ered instead of [—1,1]. The numerical method is easier to describe for the 
boundary turning point problem, thus this case also serves the purpose of 
simplifying the presentation in section 3. Subsection 3.1 introduces the dis¬ 
cretization scheme and analyzes its stability. The special mesh is described 
and the main result is stated and proved in subsection 3.2. The necessary 
changes for the interior turning point case are explained in section 4. Some 
additional remarks are also given there. This is followed by numerical results 
in section 5. 

To illustrate our main result, let us state it for the interior turning point 
case when the mesh consists of five equidistant parts (the central one around 
x = 0 being the finest). Then, an error estimate of the form 


My/I 


In 



N ~ 2 


can be proved in a discrete L 1 norm. This is an improvement over (2). 

Let us finally mention that [9] is to our knowledge the only other paper 
which uses a piecewise equidistant mesh to solve a turning point problem 
numerically. However, the problem considered there is different from (1) 
and requires a different, more complicated mesh. In particular, the number 
of equidistant parts of the mesh depends on N, which is not the case here. 


2 The Continuous Problem 

For simplicity, we are going to use e 2 instead of e in the rest of the paper. 
We consider the problem 

e 2 u" + a(x, u)v! = 0, x € I = [v, 1], u(v) = U-, it(l) = U + , (3) 

where e is a perturbation parameter, 0 < e << 1, and where v stands for 
either 0 or —1. We assume that 

a(x,u) = xb(u) (4) 
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with b £ C(K), and that U± are two different constants (otherwise U + = U- 
solves (3)). The case u = 0 describes a boundary turning point problem, 
whereas if v = — 1, we have an interior turning point problem. 

We can assume without loss of generality that £/_ < U+. Then U- and 
U + are respectively the lower and upper solutions of (3), and therefore, the 
problem has a unique solution, u e £ C 2 (I), satisfying 

u £ (x) £ U:= [U-,U+\, xel 

(see Lorenz [25]). Note that U± are also the solutions of the reduced equation 

a(x, u)u' = 0, x € /, 

subject to only one of the original boundary conditions. 

Throughout the paper, we shall assume that b £ C 3 (C/), so that u e € 
C 5 (I). Another assumption that will be needed in this paper is 

b* > b{u) > 6* > 0, u £ U, (5) 

(of course, the upper bound on b is not a restriction here). Then the solution 
u £ of the problem (3)-(5) satisfies the following estimate: 

\ u i k H x )\ < Me~ k y(x) < Me~ k z(x), k = 0,1,2, 3,4, x £ I, (6) 


where 

y(x) = e~ b ** x2 l 2£2 

and 

z( x ) = e“ m|a;|/£ . 

The above constants 6** and m are positive and independent of e. 6** 
satisfies 6** < 6*, whereas m is arbitrary. As the estimate (6) is sharp, it 
clearly shows that u e has a layer of exponential type at x = 0. Moreover, 

| u e (x) — U + 1 < My(x ) < Mz(x ), x € [0,1], (7) 


and similarly 

| u e {x) — U-\ < My(x) < Mz(x ), x € [—1,0], if v = —1. 

How to prove (6) and (7), can be found in [1] and [10], cf. [11] as well. Note 
that the proof of (6) requires b £ C 3 (U). 

It is another novelty of this paper, as compared to [1], [10], and [12], 
that the condition (5) is given locally, i.e. for u £ U, and not for u £ 1R. 
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3 The Case v = 0 


3.1 The Discretization 

Let I h denote the discretization mesh with points Xi, i = 0,1,... , N, 0 = 
xq < x\ < ■ ■ ■ < xn = 1, and let hi = x t — Xi-\, i = 1, 2,... , N. The only 
assumption on I h needed here is 

hi<h i+ 1, i = 1,2,.. rf N - 1, (8) 

(the special mesh will be introduced in the next subsection). Also, let hi = 
(hi + hi- )_i)/2, i = 1,2,... , N — 1 and x i± i = (xi + x^ ti)/2. Let w h be 
an arbitrary mesh function on I h \ {0,1}, which is identified with a column 
vector in 1R V_1 , 

w h = [wi,w 2 , • •. ,w N - 1 ] t , 

where for simplicity Wi = w^. For any mesh function, we shall formally set 
tco = U- and wjy = U + . The restriction of the continuous solution u e on 
I h \ {0,1} will be denoted by u (h Let W = {w h € 1R A _1 | vo t G U, i = 
1,2,..,, iV 1}. 

The following quantity is used to define the discretization: 

_ b*Xi-ihi 

Pi ~ 2£ 2 '* 

where b* is given in (5). Let us consider the set of indices 

J = {i\ Pi <l}^{l,2,...,N-l}. 

Note that 1 G J and let n = max J. Because of (8), 

Pi < 1, i = l,2,...,n. (9) 

If 1 < n < N — 1, we define 

(hi if 1 < i < n — 1, 

X* = \ ^ + hi+i if i = n, 

[ hi + 1 if n + l<i<N— 1. 

If n = 1, Xi = hj+i and if n = N — 1, = hi, in both cases for all 

* = 1,2,...,1V- 1. 
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We can now introduce the finite-difference operators 


D" Wi = - 
Xi 


1 f Wi- 1 - Wi w i+ 1 - w. 


h 


D'wj = 


D'wi = 


+ 

Hi 

1 - Wj~ i 

2Tn 

W’i+1 - Wi 


i +1 


= 


hi +1 

2w i+ i - Wi- Wi -1 


= 


2Xi 


The differential equation of problem (3) is discretized in the following form: 

-e 2 u" - f(x,u)' + f x (x,u) = 0, (10) 


where 


ru 

f (x, u) = / a(x,t)dt. 
Ju+ 


When discretizing / on I h , we use the notation /,; = f(xi,Wi). The notation 
fx,i, <H, bi, etc. has an analogous meaning. Then the following schemes are 
used to discretize (10): 


T c Wi = -e 2 D"wi - D'fi + f Xji , 

T +Wi = -£ 2 D" Wi - D' + fi + D°f Xti , 

T t Wi = -£ 2 D”wi - D' t fi + / Xji . 

By combining those schemes, we obtain the discretization of (10) on I h , 

Tw h = 0, (11) 


where, if 1 < n < N — 1, 

( T c Wi if 1 < i < n — 1, 

Twi = < T t Wi if i = n, 

[ T + Wi if n + l<i<N — 1. 

T c is the standard central scheme, T + is the midpoint scheme used to dis¬ 
cretize (10) at the point x i+ i, and T t is a transition scheme between T c 
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and T + . If n = 1, T c is not used and there is no need for the transition 
scheme. In that case, T = T + . Likewise, if n = N — 1, T + is not used 
and we set T = T c . For simplicity, in what follows, we shall only consider 
1 < n < N - 1. 

The central and midpoint schemes are combined above in the same way 
as in [21], but the transition scheme was not required for the type of problems 
considered there. The present transition scheme is a little simpler than the 
one used in [24]. Our stability analysis needs such a transition, but this may 
be just a technical requirement. Note that = hi+i when T + is used. This 
gives a nonstandard scheme D" for discretizing u". However, such schemes 
have been used earlier, see [12] and [20]. 

Let us finally introduce some vector and matrix norms. By || • Hoq and 
|| • ||i we denote the vector norms 

N—l 

Halloo = max \w h l Hit/bli = V IwA, 
l<i<N—1 ' 

1=1 

and, at the same time, their subordinate matrix norms. The diagonal matrix 

H = diag(xi,X 2 ,---,XiV-i) 

is used to define the following discrete L 1 norm: 

\\w h \\ H = \\Hw h \\i, 
and its subordinate matrix norm 

\\A\\ h = \\HAH~ l \\ u 

where A is an arbitrary (N — 1) X (N — 1) matrix. The vector norm 


N—l 

^2 K\wi\ 

i =1 

and the corresponding matrix norm are usually used for nonequidistant dis¬ 
cretizations of quasilinear problems, see [1] for instance. For the modified 
norms like || • \\h above, cf. [12] and [24]. 

Let G = [ gij\ = T'{w h ) be the Frechet derivative of the discrete operator 
T on mesh I h at some w h € W h . 

Lemma 1 Let (5) and (8) hold. Then G is an L-matrix. 



Proof. Since G is a tridiagonal matrix, we have to show that 


9u > 0 and g i4 ±i < 0. 

It is easy to see that ga > 0 for all the schemes used. When T c or T t are 
applied, (9) guarantees that gi.i±\ < 0. For T + , gig-\ < 0 is immediate, and 


9i,i +1 


_f!_ _ *±i + hti < bt+1 (l 

h i+ ihi h i+ 1 2 \2 


Xj+l \ 

h l+ \) 


< 0 , 


since hi + \ < Xi + \. 


□ 


Theorem 1 Let (5) and (8) hold. Then the discrete problem (11) has a 
unique solution, w^, which belongs to W. Moreover, the following stability 
inequality holds for any two mesh functions w h , v h € W: 

\\w h - v h \\ H < j-\\Tw h - Tv h \\ H . (12) 

Proof. Let e h = [1,1,..., 1 ] T € 1R A _1 . It can be shown that 

s h := e h > je h , (13) 

where the inequality should be understood componentwise. The proof of 
(13) is elementary and it requires the transition scheme. In fact, if i n + 1, 

_ K.i —1 . i K.i~h 1 ■— 7 v j 

S / — 9i—l,i ”1“ 9ii 3 9i+l,i — j 

Xi Xi 

where we formally set goi = 9n,n~ i = 0. If i = n + 1, then from T t we get 

Xi -1 „ _ £ 2 

9i—l,i — , , , i 

Xi hth l+ \ hi -|-i 


and T + gives 


and 


so that 
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The discussion above also illustrates how the transition scheme enables the 
proof of (13). 

The inequality (13) implies that G is an inverse-monotone matrix (and 
therefore an M-matrix) and also that 

IIG-'llffSy. 

This result can be applied immediately to the matrix 

A = f 1 T'(v h + s(w h — v h ))ds 

Jo 

in 

Tw h — Tv h = A(w h -v h ), 

and (12) follows. 

That (11) has a solution in W can be proved by showing that 
T(U-e h ) < 0 < T(U + e h ). 

The second inequality is immediate because of the way / is defined. To 
illustrate the proof of the first inequality, let us consider 

\T(U-e h )] n = [-Tt(U.e h )\ n = j"~ b(t) (l - dt = 

and 




(recall that wn is replaced by U + ). 

The solution is unique because of (12). □ 

3.2 The Main Result 

We shall now define the special discretization mesh. Let A be either 1/e or 
N and let 


[T+(R_e h )] JV _ 1 

2 U+ - U- f u ~ ,, , 

2 + + / m 


—e 


h 2 N 


>u+ 


XN-l 1 

hj\r 2 


< 0 


In 0 A = A, ln fc A = ln(ln fc_1 A), k = 1, 2,..., K, 
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where K = I\( A) is a positive integer such that 0 < ln A A < 1. Let £, 
1 < £ < K, denote a fixed integer independent of e. Also, let a be a positive 
constant independent of e and N. Then we define the transition points 

Tfc = ae h/ _fe+1 A, k = 1,2,... ,£. 


We shall assume that rg < 1, since e << 1. By formally setting tq = 0 and 
Tg + i = 1, we can split up the interval [0,1] into £ + 1 subintervals, 

e +1 

[0) 1] ^ k ■ Ik Tk 1 ■ T'] • 

k =1 


Each interval I k is then divided into N k > 2 equidistant subintervals, so 
that 

Ni + N 2 + ... + Ng +1 = N 


and 


N<MN k , k = 1,2,..:&£ + !. 


The points obtained in this way form the mesh on [0,1], which we shall refer 
to as the S(£) mesh (Shishkin mesh with £ transition points). The standard 
S mesh is S(l) with A = N. S(£) satisfies (8). 

Note that if (9) holds with n = N — 1 on the S(£) mesh, this practically 
means that 1/N < Me 2 , whereas usually e < 1/N. Thus it is not realistic 
to expect that the discrete operator T be identical to the central scheme T c . 
On the other hand, it is possible that n = 1 and T = T + if e is sufficiently 
small, e 2 < MN~ 2 (ln ( A) 2 . Even though we are only showing details for 
1 < n < N — 1, all our results are true for n = 1 or n = N — 1 as well. 

Let 


S = e 



2 


and 


v(x) = e~P x/£ 


(14) 


with a positive constant (5. 

We first consider the case A = 1/e. 


Lemma 2 On the S(£) mesh with A = 1/e it holds that 

^ Vi-i < M5, (15) 

where * = 1,2,... ,7V—1 and v t = v(xi) with a sufficiently large constant f3 
independent of e and N. 


11 



Proof. We shall consider several cases. If Xi € (0, t\), then h{ + 1 < 
MelV -1 h/ A and (15) follows immediately. 

In all other cases (15) can be proved with elW 2 instead of 6. If x t € 
(rfe, Tk~ i-i) for 1 < k < £, then 


h i +1 


Vi-i < 


Me 


'ln £ ~ k \" 


x * y 

MeAT 2 (ln*“ fc 


2 

e -a/31n l ~ k+1 A 
A )2-a/3 < MeN~ 2 , 


where /3 is chosen so that a/3 > 2. 

If Xi = Tki k = 1, 2,..., £, then 

Xi -‘ a Tk i 1 - V*) - T 

and (15) follows like in the previous case but with a/3 >4. □ 

Let us introduce some more notation. By X we denote 

i 

1= [jlk, 

k =0 


where 

-^k | €E {Tki Lfc+1)}• 


Lemma 3 On the S(t) mesh with A = 1/e it holds that 


iex k 


i +1 




v(x)dx < M5, 


’Xi-! 


where k = 0,1,... and where v(x) is given in ( 14 ) with a sufficiently large 
constant f3 independent of e and N. 


Proof. We have 

s m (w) 

< M ev(r k )<M5. 

The last inequality follows similarly to the proof of Lemma 2. □ 

We can now prove the main result of the paper. 
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Theorem 2 Let (5) hold and let w% be the solution of the discrete problem 
(11) on the S(i) mesh with A = 1/e. Then the following error estimate 
holds: 


— Ue \\h < MS. 


Proof. Using (12) with w h = w% and v h = uf , we see that in order to 
prove the theorem, it suffices to show that the consistency error r h = Tu% 
satisfies 

\\r h \\ H < MS. (16) 

Let 

r % = T c v!) and = T + v!f. 

Since the transition between T c and T + can occur at any mesh point, we 
can consider separately ||r^||// and ||r+||#. As for T), it is only used at x n 
and therefore, it is sufficient to estimate x n |r n |. 

When the central scheme is used, Taylor’s expansion of r t = r C} i about 
Xi gives 


Xi\n\ < M \ e 2 


^L+i(^L+i hi)\v^. \xi)\ + h 


i +1 


r%i +1 


'Xi -1 


u^\x)\dx 


rXi- |_i | 

+ h i+ i(h i+ i - hi)\p"(xi)\ + hf +1 I \p (3) (x)\dx \ , 

Xi-1 


where p(x) = f(x,u e (x )) (for the integral form of the consistency error, cf. 
Kellogg and Tsan [26] for instance). Then (6) implies 


Xi\n\ < M(Qi + Ri ), 


where 

and 


Qi — ^h+i(^h+i hi)e Vi—i 

(^r 1 ) LZ‘ v(x)dx - 


Ri = 


Here, v(x) is like in (14), with a constant (3 satisfying 0 < (3 < m. (3 
is independent of e and N, it can be chosen arbitrarily close to m, and 
therefore it can be made sufficiently large. On writing 


r h \\H = XiN + ^2xi\ri\, 

igz ieX 
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from Lemma 3 we immediately get 


Y,Xi\ri\<MY,Ri<M6. 
iei iei 

On the other hand, if i 0 X, i.e. if xi = 77 it holds that 

h? 

Qi < -^Vi -1 < MS, 

£ 

because of Lemma 2 . Since the number t of transition points does not 
depend on N, this proves (16) for r h = r^. 

Let us now consider the consistency error when T = T + . Then, r, = r +) j 
is expanded about x i+1 / 2 - The errors due to D' + and D° can be treated as 
above, but the one due to D" requires a closer attention. Let 

r” = £ 2 [D"u e (xi) - u"(x i+ 1 / 2 )]. 


We have 


XiWl | < Me 2 


r^i -\-1 

(^■i+1 (•^ 1 + 1 / 2 ) I T hi +1 / 

Jxi—1 


,(3) 


(x)|dx 


Since in this case pi > 1, it holds that 

e 2 < Mxi-ihi. 


This and ( 6 ) imply 


Qi — hi(h %+1 hQe Vi— 1 , 
i?" = /l 2 +i ^zl [ X ‘ +1 z ( x ) dx . 


•Xi~ 1 


We can handle Q" in the same way as Qi above, and 

'/?. i+ i\ 2 Xj_i _ 


K < 


< M 


i +1 


mxi-i/2e f 4+1 e -mx/2e^ x 
Jxi-i 


r +1 e~ mx/2£ dx , 

Jxi-i 


so that Lemma 3 can be applied. Thus, (16) holds with = r 


h 
+ * 
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Finally, let as consider the transition scheme at x n . We again expand r n 
about x n to get 


Xn\r n | < M<s 


h n +i\u"(x n )\ + h 2 n+1 max |4 3) (^)l 

Xi-l<X<Xi + l 


(17) 


+ h-n+ilp" (0)\} 

with some 6 € (x n _i, x n+ \). The fact that in this case p n +1 > 1 implies 


s £ Mx n h n -\. i* 


We need this inequality and (6) in the estimate 

h n +i\u"(x n )\ < Mhl +1 ^z(x n ) < M —■ v n . 

The other terms on the right-hand side of (17) can be estimated directly 
using (6). Thus, it follows that 

h 2 

Xn\r n \ < M-^Vn-l < MS , 

where the last inequality follows from Lemma 2. □ 

Let us now turn to the case A = N. 


Theorem 3 Let (5) hold and let w% be the solution of the discrete problem 
(11) on the S(£) mesh with A = N. Then the following error estimate holds: 

|| w h e - u h £ \\ H < M (5 + lV-( 1+??lnJV )) , 

where q is some positive constant independent of e and N. 


Proof. The result can be proved analogously to the proof of Theorem 
2, cf. [18] as well. The only case which produces the lV^( 1+r?lnjV ' ) term is 
when Xi = T£, T£ + k, where n is the mesh step in [re, 1]. For instance, when 
estimating expressions like the left-hand side of (15) (which is also needed 
in the proof of Lemma 3), the best we can get at those points is 


h 2 


■j+i 


Vi -1 < M——e~ a ^ nN = M 
N Z £ 


1 


eN 2 + a P' 


15 



The case x t G (rg + k, 1) still gives the same result as before, but requires a 
somewhat different technique: 


h i +1 ^ £ 

T -Vi -1 < M 


1 


£ 


TV 2 N a Pe 2 


; -7/iV£ < M 


TV 2 ’ 


where 7 is a positive constant independent of e and TV, and where a(3 > 2. 
Therefore, when x t = rg, Tg + K , we use a different estimate of the consistency 
error: 


Xi\n\ = Xi\Tu e (xi)\ < Mh i+ i{e 2 max [\u''(x)\ + \p'(x)\ + |g(s)|]}, 

Xi-l<X<Xi + l 

where q(x ) = f x (x,u e (x)). We use here estimates ( 6 ) and (7) with y(x) to 
get 

Xi\n\ < Mh i+ iy(xi-i) 

and 

X i\n\ < MN~ 1 e~‘ n ^ nN ^ 2 < MN~( 1+rtlnN \ 

on setting i] = b^a 2 / 2 . □ 

The error estimate of Theorem 2 is better than that of Theorem 3 when 
e —>• 0. It may look like the term TV^ 1+,?ln in the estimate above arises 
for purely technical reasons. However, the numerical experiments in section 
5 will show that when that term dominates in the error, the error does not 
decrease together with e. 


4 The Case v = — 1 and Other Remarks 


Theorem 2 holds also for the case v = — 1, when the problem (3) is considered 
on I = [— 1 , 1 ], if the numerical method is modified appropriately. As for 
the mesh, the simplest thing to do is to extend S(£) to [— 1 , 0] symmetrically 
to xo = 0: 

x-i = -Xi, i = 1,2,... ,TV, 
h-i = x-i - x-i-i = h i+ 1 , i = 0,1,. .., TV - 1. 

This is accompanied with other symmetrically changed definitions, like 


Xi = < 


hi 

u I hi +1 
rij H-7j— 

hi 

+ /q +1 
^■i+l 


if — TV + 1 < i < —n — 1, 
if i = — n, 

if — n + 1 < i < n — 1 , 
if i = n, 

if n + l<i<N — 1, 
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and 


D'_Wj = 


Wj - Wj-l 
hi 


D[_Wi = 


Wi+i +Wi - 2wi-i 


D°_ Wi = 

The scheme is also symmetrical: 


2 Xi 

Wi + Wi-l 


Twi = { 


T-Wi 

T t -Wi 

T c Wi 

T t Wi 

T+Wi 


if — N + 1 < i < —n — 1 
if i = — n, 

if — n + l<i<n — 1, 
if i = n, 

if n + l<i<N — 1, 


where 

T-Wi = -£ 2 D" Wi - D’_fi + D°_fx,i, 

and 

T t -Wi = - e 2 D"wi - D[_fi + / X)i . 

Note that T c is always used at least at xo = 0. To make the discretization 
even more symmetric, we also replace 17+ in / with U- when discretizing 
(10) at Xi € (—1, 0), and with (17+ + 17_)/2 at xo = 0. 

Under the conditions of Theorem 2, the following estimate holds: 


\w h £ -u h E \ 


oo — 


M 


h/(l/e) 

N 


(18) 


This is because the smallest mesh step is hi = t\/N. The above result does 
not mean e-uniform pointwise convergence, but hr (1/e) grows very slowly 
when e —> 0. For instance, if i = 3 and e = 10 -12 , h/(l/e) ~ 1.2. The 
numerical results of the next section show that pointwise accuracy is even 
better than what (18) indicates. A similar discussion is true for S(£) with 
A = N. 

It is possible to adjust the present method to the most general problem 
(1) under the conditions given in [1], The reduced solutions are not 17_ and 
17+ in that case, but some more complicated functions. They have to be 
incorporated in the function /, like it was done in [1], In order to prove the 
stability inequality corresponding to (12), the conditions of type (5) on b 
have to hold for u € 1R. The same has to be assumed of all conditions on 
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c. This is because in this case it is generally difficult to find the upper and 
lower solutions of (1) and of the discrete problem. For Theorems 2 and (3), 
more complicated estimates corresponding to ( 6 ) and (7) have to be used, 
see [ 1 ]. 


5 Numerical Results 

Our test problem is more general than (3) but less general than (1), 

—e 2 u" — xuu' + c{x) = 0, tt(zbl) = U±, 

where c(x), U- ~ 1, and U + ~ 3 are determined by the exact solution being 

x 

u e (x) = 2 + tanh —. 

We have tested the S(£) mesh with both A = 1/e and A = N, and with 
values of £ = 1, 2, 3. Let 

N k 

Qk = ~jy, k = 1,2, ...,£+ 1, 

where N is the number of mesh steps in [0,1]. Table 1 shows the values of 
the ratios qk that are used in the meshes below. They were kept fixed for 
each £ regardless of other mesh parameters, including N. 


Table 1. Mesh parameters 


£ 

1 

2 

3 

qi 

3/4 

1/4 

1/8 

Q2 

- 

1/2 

1/8 

73 



1/2 


The choice of the ratios qk may influence the errors significantly. Other 
ratios have been also tested and the results for this test problem were the 
best when there were around 75% of the mesh points in the layer. The 
question of the optimal choice of the ratios is problem-dependent and seems 
to be difficult to solve in general. 

The tables below show the errors 

E = E(N) = \\ W £ — u^Hoo and Ex = E^N) = \\w% - u h e \\ H , 


18 




where the dependence on N indicates that the mesh is used with 2N mesh 
steps in [—1,1]. The corresponding numerical orders of convergence, Ord= 
Ord(lV) and Ordi=Ordi(lV), are also listed, where 


Ord (AO = 

and Ordi(A^) is defined analogously. The results in Tables 2-4 are given for 
the transition point coefficient a = 1, whereas Table 5 presents some results 
for a = 2. 

Table 2 illustrates that the errors are smaller if £ is larger, that is, if the 
meshes have more subintervals within the layer. This is to be expected. The 
case A = 1/e is shown. If A = N, the errors behave analogously with respect 
to the change in l. 


Table 2. S(£) mesh with A = 1/e, a = 1, N = 512 


e 

£ = 


£ = 

2 =ts 

£ = 

2=21 


i 

7.94-05 

1.19-08 

1.31-04 

1.06-09 

1.96-04 

9.83-11 

E 

Ei 


2.00 

1.89 

2.00 

1.99 

2.00 

2.00 

Ord 

Ordi 

2 

3.84-05 

6.91-09 

4.74-05 

4.15-10 

5.53-05 

3.00-11 




2.00 

1.81 

2.00 

1.98 

2.00 

2.00 



3 

1.96-05 

4.29-09 

2.29-05 

1.93-10 

2.60-05 

1.33-11 




2.00 

1.68 

2.00 

1.96 

2.00 

2.00 




In the remaining tables, only £ = 3 is considered. Table 3 gives more 
details of what can already be observed in Table 2, viz. errors E\ decreasing 
together with e, while at the same time, errors E slightly increase, still 
preserving a high accuracy. The increase of E is what we can expect from 
(18), but the accuracy and its order are higher than what (18) indicates. 
Both Ord and Ordi are around 2 , with Ordi being somewhat smaller for 
larger e values. 
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Table 3. S(3) mesh with A = 1/e, a = 1 


N 

£ = 

2= T1 

£ = 

2 = tb 

£ = 

2^2 


64 

1.18-03 

1.88-07 

1.57-03 

1.73-08 

1.73-03 

1.27-09 

E 

Ei 


2.01 

2.06 

1.51 

2.00 

1.48 

1.89 

Ord 

Ordi 

128 

3.13-04 

4.80-08 

3.63-04 

2.97-09 

3.94-04 

2.45-10 




1.92 

1.97 

2.12 

2.54 

2.13 

2.37 



256 

7.83-05 

1.38-08 

9.16-05 

7.51 10 

1.04-04 

5.33-11 




2.00 

1.80 

1.99 

1.98 

1.92 

2.20 



512 

1.96-05 

4.29-09 

2.29-05 

1.93-10 

2.60-05 

1.33-11 




2.00 

1.68 

2.00 

1.96 

2.00 

2.00 




When A = N, both E and E± errors are more uniform in e but the 
accuracy is worse than for A = 1/e. This is shown in Table 4, where we 
can see that Ord is less than 2 as can be expected from the error estimate 
of Theorem 3. Also, E\ is much worse than in Table 3 even though Table 
4 shows Ordi significantly higher than 2. In fact, in Table 4, E\ does not 
decrease together with e, which is possible according to the error estimate 
of Theorem 3. This can be improved by increasing the value of a so that 
the ^-independent, term of the error estimate gets negligible. Table 5 shows 
that E\ errors decrease for a = 2 when £->0 and that they become almost 
as accurate as in Table 3, but at the same time, the larger a spoils E a little. 


Table 4. S(3) mesh with A = N, a = 1 


N 

£ = 

2= T3 

e = 2- 

18 ^ n 


64 

1.89-03 

9.21-05 

1.90-03 

9.21-05 

E 

E\ 


1.91 

2.62 

1.91 

2.62 

Ord 

Ordi 

128 

5.07-04 

1.41-05 

5.10-04 

1.42-05 




1.90 

2.70 

1.89 

2.70 



256 

1.37-04 

2.10-06 

1.39-04 

2.10-06 




1.89 

2.76 

1.88 

2.75 



512 

3.64-05 

3.02-07 

3.76-05 

3.04-07 




1.91 

2.79 

1.88 

2.79 
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Table 5. S(3) mesh with A = N, a = 2 


N 

£ = 

2= T1 

£ = 

2 = tb 

£ = 

2T 72 


64 

1.59-03 

2.63-07 

1.59-03 

3.02-08 

1.59-03 

1.56-08 

E 

Ei 


3.09 

3.15 

.09 

4.07 

3.09 

4.66 

Ord 

Ordi 

128 

4.44-04 

4.81-08 

4.44-04 

3.51-09 

4.44-04 

7.24-10 




1.84 

2.45 

1.84 

3.10 

1.84 

4.43 



256 

1.40-04 

1.67-08 

1.40-04 

1.07-09 

1.40-04 

8.44-11 




1.66 

1.53 

1.66 

1.73 

1.66 

3.10 



512 

4.34-05 

5.53-09 

4.34-05 

3.46-10 

4.34-05 

2.23-11 




1.69 

1.59 

1.69 

1.61 

1.69 

1.92 




How fast do the E\ errors in Tables 3 and 5 decrease as £ -A 0? This 
e-order can also be measured numerically, analogously to Ordi: 

£,(40 - Ei(s) 

U1CI '“ In 4 

where N is kept fixed and the dependence of E\ on the value of £ is expressed. 
The results are given in Table 6 . We can see that they confirm the expected 
value of 1 for A = N and a = 2, whereas for A = l/£, they are better than 
what Theorem 2 indicates, particularly for larger e. 


Table 6. Ord e on S(3) mesh with N = 512 


£ 

A = 1/e, a = 1 

A = N, a = 2 

2-u 

1.45 

1.00 

2~ w 

1.21 

1.00 

2 -ls 

1.03 

1.00 

2 ” 2U 

0.97 

1.00 

2 -22 

0.96 

0.98 


Let us finally mention that we have also tested the discretization in which 
the central scheme is used instead of the transition schemes Tt and T)_. The 
errors are somewhat worse, the difference being greater when A = 1/e than 
when A = N. Thus, it may be possible that the use of the transition schemes 
is not entirely for technical reasons. 
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